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Coordinate transformations are powerful tools for the solution of the 
partial differential equations which describe physical phenomena. The 
use of transformations leads to well ordered discretizations of the 
physical domain and thereby renders a simplification in a numerical solution 
process. The discretization is constrained by the underlying physics, the 
problem geometry and the topology of the region where the solution is to 
be obtained. The constraints can be stated in geometric terms. In 
particular they can be categorized as boundary constraints, uniformity 
constraints, and internal constraints. Boundary constraints include: the 
basic geometry of solid objects, the transmissive junctures between and 
around solid objects, the choice of representation for the boundaries, the 
angles at which transverse coordinate curves intersect boundaries, and the 
rate of entry for such coordinate curves. Uniformity constraints are 
applied to either local or global distributions of coordinate curves or 
points to form a basis from which the curves or points can be redistributed. 
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This may be based on an a priori specification of a distribution function 
or on a solution adaptive approach. In either case, the redistribution 
must not be distorted by the underlying transformation. Internal con- 
straints are applicable when an interior shape or interior mesh structure 
is to be smoothly embedded within a global mesh to simplify the simulation 
of physical processes in the given region. 

Algebraic mesh generation techniques are highly advantageous for 
meeting the constraints described above. Algebraic techniques provide 
exact control of the mesh properties necessary to satisfy the given 
constraints. Although other methods have been developed which provide some 
degree of control, the level of control is not in general sufficient to 
satisfy certain of the constraints. For example, the smooth embedding of 
a Cartesian mesh within a global mesh structure cannot be readily con- 
structed with the application of differential equation techniques. Also, 
three dimensional meshes are not in general readily obtained with non- 
algebraic techniques. On the other hand, algebraic techniques require 
more complex specification of the data to assemble a mesh. The purpose 
of this paper is to present an overview of algebraic techniques for mesh 
generation and set forth the underlying concepts which have been successful. 
Both two- and three-dimensional domains are considered. 

The Multi-Surface Transformati on 

When curvilinear coordinates are employed in the numerical solution 
of a boundary value problem, constraints must often be placed upon the 
coordinates, in addition to the basic requirement that the bounding sur- 
faces are coordinate surfaces of one or more coordinate systems. The 



locations of the constraints can occur anywhere in the problem domain. 

On the boundaries, a particular pointwise distribution may be needed; in 
regions near boundaries, a particular coordinate form may be advantageous; 
and away from the boundaries, an internal coordinate specification may 
also be required. Typically, the constraints will arise either to resolve 
regions with large solution gradients or to cause some simplification in 
the problem formulation and solution. 

In conjunction with the demand for constraints, the general multi - 
surface transformation [1] will be examined. The multi-surface transfor- 
mation is a method for coordinate generation between an inner bounding 
surface ^ and outer bounding surface P„. To establish a particular 
distribution of mesh points on each bounding surface, a common parameteriza- 
tion t is chosen for each surface. This is equivalent to a coordinate 
description of the surfaces which yields the desired surface mesh when the 
parametric components of t are given a uni form discretization. With the 
parametric description, the inner and outer bounding surfaces are denoted 
by and ? N (t) respectively. In continuation, parameterized 

intermediate surfaces f^(t) » • • . (t) are introduced so that they 

can be used as controls over the internal form of the coordinates. The 
intermediate surfaces are not coordinate surfaces but, instead, are 
surfaces which are used to establish a vector field that is composed of 
tangent vectors to the coordinate curves spanning the coordinate system 
to connect bounding surfaces. It is also assumed that the collection 
of surfaces F^(t), » • • • »^(^) is ordered from bounding surface 

to bounding surface. An illustration is given in Fig. 1. For a fixed 
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Fig. 1 - A piecewise linear curve and its tangent field. 


parameter value t, there is a corresponding point on each surface. The 
piecewise linear curve obtained by connecting corresponding points is 
given by the dashed curve in Fig. 1. From the figure, it can be observed 
that the tangent directions determined by the piecewise linear curve are 
piecewise constants. As t is varied, the field of tangent directions 
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w' obtain their smoothness (level of differentiability) only in t. To obtain 
smoothness in going from bounding surface to bounding surface, a suffi- 
ciently smooth interpolation must be performed. The result is a smooth 
vector field of undetermined magnitude which gives the desired tangential 
directions for coordinate curves connecting the bounding surfaces. A 
unique vector field of tangents is then obtained by correctly choosing 
magnitudes so that, on integration, the bounding surfaces are fit precisely. 

In symbols, a vector field tangent to the piecewise linear curves is 
given by 


V*> \^k+l^ ' ’ 


( 1 ) 


between the kth and (k + 1 ) th surfaces where k is taken to vary (if 
N > 2) from the first bounding surface to the final intermediate surface. 
These vectors are indicated in Fig. 1. The coefficients are scalars 
which determine the magnitude of the vectors but not the directions. When 
an independent variable r is assumed for the spanning direction, a 
partition r-| < ... < r^_-j can be specified in correspondence with the 
tangents of Eq. 1. The partitioned variable can then be used to represent 
the tangents as the discrete vector valued function which maps r^ into 
for k = 1 ... N-l. A sufficiently smooth vector field tf(r,t) is 
then obtained by a sufficiently smooth interpolation ^(r^t) = ^(t). 

With r as a continuous independent variable, the r-derivative of the 
coordinate transformation ^(r,^) is equal to the interpolant and is 
given by 
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( 2 ) 


I = f ■ £ \< r A& ■ 

where i\(r.) is unity at k = j and vanishes otherwise. When Eq. 2 

K J 

is integrated with an initial r-| value of ^(t), the transformation 
becomes 




^(r,t) = ^(t) + 


N-l 

Z 

k=l 


w^Vi^ * 


( 3 ) 


where 


G k (r) = / ^ k (x)dx , (4) 

r l 

from which we observe that the interpolants il> k must be continuous!” 
differentiable up to an order which is one less than the level of smoothness 
desired for the coordinates. The construction of local controls on the 
coordinates will rely upon the development of suitably smooth interpolation 
functions. If the magnitudes A^ of Eq. 1 are chosen so that each 
A k G k (r N _i) is unity, then the evaluation of the transformation at r^ -j 
will reduce to f^(t) by means of a telescopic collapse of terms in 
Eq. 3. With this choice, we obtain the general multisurface transformation 
of Eiseman [1] which is given by 

N-l G.(r) 

?(r,t) = ?,(t) + I Trf- T [i* k+ 1 (T) - ? k (T)] . ( 5 ) 

1 k=l b k lr N-l ; K 1 K 


V* 
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On examination, each interpolation function ^ can be rescaled without 
changing the transformation; hence, the original vector field interpolation 
becomes an interpolation only on vector directions. 

When the interpolants ^ are polynomials in r, the coordinate 
curves which connect the bounding surfaces are globally defined by poly- 
nomials in r of one greater degree. The specification of boundary 
properties for the curves and a global control over their general form are 
obtained by choices of intermediate surfaces and the associated partitions 
of r. For notational simplicity, let r-j = 0 and r^_-j =1. In the 
simplest case when there are no intermediate control surfaces, there is 
just one vector field direction ^(t) which is determined solely by 
the bounding surfaces. The interpolant tjj-j is then a constant function, 
G-j(r) = rip-j, G-j(r)/Gj(l) = r, and the polynomial 2-surface transformation 
becomes 

?(r,l) = fyt) + r[? 2 (t) - ?,(*)] , (6) 

which is the case of linear coordinate curves connecting boundaries. The 
linear case has occurred in many studies including [2], [3], and [4] and 
is limited to at most one prescribed coordinate property per boundary 
which can be either a pointwise distribution or a distribution of angles 
with the linear transverse coordinate curves. To allow for the prescrip- 
tion of an additional coordinate property on one of the boundaries, an 
intermediate control surface is introduced and the polynomial 3-surface 



transformation is computed from Eq. 3 with if>i * 1 - r and = r 
corresponding to directions of \^(t) and Eq. 1 and Fig. 1. 

The integrals from Eq. 4 become 




G ] (r) = r - ~ , 

(7) 

2 

S 2 (r) = — , 


and since G^(l) = G^O) = the original vector field which is discrete 
in r is determined by 


f,(t) - 2[P 2 (t) - P,(t)], 
? 2 (t) - 2[? 3 (t) - ? 2 (t)]. 


( 8 ) 


because = l./G^O) in Eq. 1. Upon substitution from Eq. 7 the 
polynomial 3-surface transformation is given by 


P(r.l) • f,(t) + r(2 - r)[P 2 (t) - P,(t)] 
+ r 2 [f 3 (t) - P 2 (t)]. 


(9) 


In continuation, an additional coordinate property can be prescribed 
on each boundary when two intermediate control surfaces are used. The 
polynomial 4-surface transformation is computed from Eq. 3 with interpolants 

V 
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^(r) = (1 - r)(r 2 - r) , 


^ 2 (r) = r( 1 - r), (10) 

i^ 3 (r) = (r - r 2 )r, 

which respectively correspond to the directions of \f-|(t), \f 2 (t), and 

^ 3 (t) which, in turn, are respectively associated with partition points 
r-| = 0, r 2 , and r 3 = 1 and which are defined to vanish at all partition 

points except the ones of association for each function. For simplicity, 
we will set r 2 = j so that the partition is uniform. The nonuniform 
case is a simple but algebraically more complex extension [1]. With 
r 2 = ■£-, the integrals from Eq. 4 become 

Gi(r)=lr-fr 2 + lr 3 , 

G 2 (r) = \ r 2 - 1 r 3 , (11) 

G 3 (r) 4r 3 -}r 2 , 

and from an evaluation at the endpoint r = 1 we have G^ ( 1 ) = jy, 

G 2 (l) = and G 3 (l) = By substitution, the polynomial 4-surface 
transformation is given by 
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£(r,T) = ^(t) + r(6 - 9r + 4r 2 )[? 2 (t) - ^(t)] 
+ r 2 (3 - 2r)[? 3 ('{) - ? 2 (t)] 

+ r 2 (4r - 3)[? 4 (t) - ? 3 (t)], 

with r-derivative 

|” = 6(1 - r) ( 1 - 2r)|? 2 (t) - ^(t)] 

+ 6r(l - r)[? 3 (t) - fyfy] 

+ 6r(2r - l)[? 4 (t) - ? 3 (t)] . 

By direct evaluation 

?(0,t) = 

?(l,t) = fyt), 

and 

|E (0,t) = 6[? 2 (t) - ?,(*)]. 

J~ (l,t) = 6tfyt) - ? 3 (t)]. 




( 12 ) 


(13) 


(14) 


05) 


% 
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which explicitly shows that in addition to fitting the boundaries (Eq. 14), 
the intermediate surfaces i^(t) and F^(t) can be used to control the 
angles at which the transverse coordinate curves in r intersect the 
boundaries (Eq. 15). Moreover, the choice of intermediate surfaces can 
also be used to control the shape of the transverse curves and the distribu- 
tion of the constant r coordinate surfaces. The general derivation and 
discussion is given in Eiseman [1], For our purposes, the discussion on 
coordinate system controls will be deferred until local methods are pre- 
sented for our survey of some of the material developed in Eiseman [5], 

[ 6 ]. 

An alternative form of the polynomial 4-surface transformation (Eq. 12) 
can be obtained from the evaluations of the transformation (Eq. 14) and its 
derivatives (Eq. 15). With the evaluations, the intermediate surfaces 
can be expressed entirely in terms of boundary data, which results in 

? 2 (t) - f(O.t) + (O.t), 

(16) 

and 

P 3 <t) = f(i.t) - O.t). 

Upon substitution of Eqs. 14 and 16 into Eq. 12 we obtain 


?(r,t) = (1 - 3r 2 + 2r 3 )£(0,t) + r 2 (3 - 2r)?(l,t) 

+ r(l - r) 2 || (O.t) - r 2 (l - r) (l,t). 


(17) 
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after grouping terms by boundary type. By examination, the coefficients 
of the boundary evaluations for the transformation and its r-derivative 
can easily be identified as just the canonical Hermite cubic interpolants 
on the unit interval 0 <_ r <_ 1 . When the r-deri vati ves are specified 
to be normal to the respective boundaries, we obtain the transformation 
presented by Smith and Weigel [3]. 

In continuation, polynomial N-surface transformations can be 
systematically established from Eqs. 3-4 and the interpolants 

N-l 

iMr) = n (r - r,), (18) 

K i = l 1 

i*k 

for k = 1,2,..., N-l. In each case, the transverse coordinate curves are 

polynomials of degree N-l in r with vector valued coefficients which are 
functions of the surface coordinates t. Polynomials, however, are globally 
defined for all r, and as a consequence, local mesh properties cannot be 
controlled without a global effect. As an example, suppose that we wish 
to smoothly embed a general rectilinear Cartesian system within a global 
mesh structure to obtain a system of the form illustrated in Fig. 2 where 
the Cartesian region within the mesh is bounded by the darkened curves. 

In the Cartesian part of the mesh, coordinate curves in r would be lines 
which pass through it at a uniform rate. Since global polynomials in r 
would be uniquely determined by the Cartesian region, curved boundaries 
could not be fitted. 
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Fig. 2 - A smoothly embedded Cartesian region 
within a global mesh structure. 

To obtain precise local controls which could be successfully applied 
to generate a mesh as illustrated above in Fig. 2, local forms of the 
multisurface transformation (Eqs.3-4) were established and analyzed by 
Eiseman [5], [6], and [7]. Our discussion will follow the development 
given by Eiseman in [6] and will focus upon two-dimensional applications 
with a surface coordinate t = t. When the interpolants 4^ are 
nonvanishing on only a local region, the precise local controls over the 
coordinates that were obtained will be illustrated with the local piecewise 
linear interpolants that are depicted in Fig. 3. For algebraic simplicity, 
the analysis is restricted to the case with the uniform partition r^ = k 
with the clear understanding that nonuniform partitions will follow the 
same analytic pattern. Since the multi-surface transformation remains 
unchanged when the interpolants are scaled by any sequence of nonzero 
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Fig. 3 - Piecewise linear local interpolants with partition points 

r^ = k for k = 1 ,2 , . . . ,N-1 . 


numbers, the height ^(r^) of each interpolant can be chosen arbitrarily. 
In particular, the form of the multi-surface transformation can be 
simplified when the heights are adjusted so that each interpolant integrates 
to unity which yields ^(r^i) = 1 for all k. The integrals are 
obtained from triangular areas, and by direct observation, lead to the 
height adjustments ^(r-j) = 2, ip k (r k ) = 1, and <I' N _ 1 (r N _ 1 ) = 2 in 
correspondence with the successive illustrations in Fig. 3. Also, in 
correspondence, the explicit form of the normalized interpolants are 
given by 


^(r) 


2(2 - r) 


for 1 < r < 2 


for 2 < r < N - 1 


(19) 


S 
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r 


O k (r) = 


0 for 1 £ r < k - 1 

(r - k) + 1 for k - 1 £ r < k 

(k - r) + 1 for k £ r < k + 1 

0 for k + 1 < r < N 




> ; ( 20 ) 


- 1 


^N_l( r ) = < 


for 1 < r < N - 2 


2(r - N + 2) for N - 2 £ r £ N - 1 


and their integrals defined in Eq. 4, by 


(21) 


Gi(r) 


1 - (2 - rr for 1 £ r < 2 


for 2 < r < N - 1 


( 22 ) 


G k (r) 


0 




l/2(r - kr + (r - k) + 1/2 
-l/2(k - r) 2 - (k - r) + 1/2 


1 


for 1 £ r < k - 1 

for k - 1 £ r < k 

for k £ r < k + 1 

for k + 1 < r < N 


>5(23) 


- 1 


Si-1 ( r ) 


for 1 < r < N - 2 


(r - N + 2) for N - 2 £ r £ N - 1 


(24) 


which are depicted in Fig. 4. 
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Fig. 4 - Integrals of normalized interpolants for the partition r k = k. 

On the interval k £ r < k + 1, the integrals G^(r) which correspond 
to interpolants defined over nonintersecting intervals are either unity or 
vanishing depending upon whether the interval of definition precedes or 
follows the interval under examination. When i = l,2,...,k - 1 which 
is nonvacuous for k > 1, the integrals G^r) have been evaluated over 
the entire domain for which the respective interpolant i|^ is nonvanishing; 
hence, these preceding integrals are unity by the chosen normalization. 

When i = k + 2, k + 3, . . . ,N - 1 , the interpolants each vanish on 
1 £ r < k + 1; hence, the integrals G^ also vanish there. As a 
consequence, G^ for i = k, k + 1 yield the only nontrivial contribu- 
tions for the multi-surface transformation which reduces to 

f(r.t) - ?,(t) + £ [? 1+] (t) - P,(t)3 + G k (r)[? k+1 (t) - ? k (t)] 

+ G k+1 ^ r ^k+2^ " ^k+1 (t)] + 0 (25) 

= \(t) + G k( r ) [^k+1 ( t ) " + G k+l^ r ^k+2^ ’ ^k+1^^’ s 

88 




which depends upon only the three control surfaces \+] » \+2 

which can be arbitrarily selected to our advantage when they are not 
bounding surfaces. With substitutions from Eqs. 22-24, we obtain the 
partition point (r^ = k for k = 1,2,...,N - 1) evaluations 


?o,t) =■ p,(t), 

P(2,t) = \ [f 2 (t) + f 3 (t)], 

P(k,t) - \ [f k (t) + f k+ i(t)]. (26) 

\ c^ N _ 2 ( t ) + 

?(N-l,t) = ? N (t), 

which, in addition to boundary fitting at the end points r = 1 and 
r = N - 1, also shows that the transformation passes through the midpoints 
of the lines which connect the intermediate control surfaces for any 
fixed surface coordinate t. Moreover, from the general multi-surface 
construction (Eqs. 1-5 and Fig. 1), the transverse coordinate curves are 
tangent to the connecting lines at the partition point evaluations. The 
tangents at partition points can be explicitly obtained from substitutions 


?(k+l,t) = 


?(N-2,t) = 
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of the interpolation functions (Eqs. 19-21) into the r-derivative of 
the transformation (Eq. 25) and are given by 


dP 

3r 

3r 


9 ? 

3r 

3 ? 

3r 


3 ? 

3r 

3 ? 

3r 

In graphical 


O.t) = 2[P 2 (t) - ^(t)]. 

(2,t) = P 3 (t) - ? 2 (t), 

(k,t) = ? k+1 (t) - ? k (t), 

{k+1 »t) = ? k+2 (t) - ? k+1 (t), 

(N-2,t) = ? N _.,(t) - P N _ 2 (t), 

(N-l,t) = 2[? N (t) - ? N _i(t:)]. 

form, this process is depicted in Fig. 5. 



I! M | li|!! If'"' 



Fig. 5 - Coordinate curve segments for k < r < k+1. 


Between the control surfaces P. and P. for i > j, the distribution 

• J 

of constant r coordinate surfaces can be controlled for the general 
multi surface transformation (Eqs. 4-5) when uniformity can be specified 
along a direction of measurement 


i„) . M! -W- 

iifyt> - 


(28) 


for then arbitrary distributions can be applied relative to uniform 
conditions. An illustration is given in Fig. 6. 
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Fig. 6 - Measurement of uniformity. 


To obtain uniformity, the projected arc length 

S p (r,t) = [?(r,t) - fyt)] • i(t), (29) 

depicted in Fig. 6 must be linear in r, or equivalently 9Sp/3r must be 
independent of r. But then from Eq. 25 and with the relative projections 
C m (t) = [? m+1 (t) - i* m (t)] • ?(t), we have 

9S 

KZp' = *l J L,(''')C|.(t) + 'i’l.j.-i (*' , )Cr, + i (t) x 



' -2C,(t) + C 2 (t) 

’ < - C k (t > + c k-l (t > 

I + ^N-l^ 

where the last equality comes from Eqs. 19-21. Hence, for k = j, j+1,..., i-1, 
uniformity is obtained if 2C-|(t) = C 2 (t) should k = 1, ( t ) = ( t ) 

should 1 < k < N - 1, and C^ft) = 2C^_-|(t) should k = N - 1. A more 
thorough discussion on uniformity is available in Eiseman [1] for the 
global case, in Eiseman [5], [6] for the local case, and in Eiseman [7] 
for the general cases. 

To explicitly demonstrate the application of the local controls, and 
at the same time, reveal the basic algorithmic steps, coordinates will 
be obtained for a simple transition from a purely Cartesian system into a 
purely Polar system. For 0 £ t <_ 1 , the Cartesian coordinates will be 
specified below a line $(t) = (2t-l,0) and the Polar coordinates, beyond a 
circular arc /2 u(t) where u(t) = (-cos o, sin e) for o = (2t + l)n/4. 

The line and the arc are depicted in Fig. 7. 


for k = 1 

for 1 < k < N - 1 / r+ function of t, (30) 
for k = N - 1 
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Fig. 7 - Basic curves for the Cartesian to Polar transition example. 


To obtain uniformity near the sides (t = 0,1) of the transitional region, 
the unit vertical distance will be used as a basis for displacements to 
establish uniformity below the line and beyond the circular arc. For the 
line, let 


^(t) = $(t) - ( 0 ,f), 

P 2 (t) = Q(t) - (0,2), 
? 3 (t) = <$(t) - (0,1), 
P 4 (t) = Q(t) , 



I ■ !'!*!• 


so that for 


T(t) « T" = (0,1), 

IlfyO - 


( 32 ) 


we have 


C,(t) - [fyt) - ?,(t)] • f(t) = (0,1) • (0,1) . 1, 

C 2 (t) = [? 3 (t) - ^2^)3 • T(t) = (0,1) • (0,1) = 1, (33) 

C 3 (t) = [? 4 {t) - ? 3 (t)] • ^t) = (0,1) • (0,1) = 1, 

which satisfies uniformity for 1 £ r <_ 3 and yields a Cartesian system from 
?(l,t) = ^(t) = §(t) - (0,|), (34) 


up to 


P(3,t) * 1 [P 3 (t) + P 4 (tn = 5(0 - (0,1). 


(35) 


Similarly, for the circular arc, let 


? 5 (t) = /2 u(t), 
fyt) « (1 + /2) u(t). 


(36) 
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? 7 (t) = (f + /2) G(t), 


be the last three surfaces so that for 


fyt) - ? 5 (t) 

ifyt) - ^Tt)Tl 


u(t) 


we have 


(37) 


C 5 (t) = [£ g (t) - j^ 5 (t)] • x(t) = u(t) ♦ u(t) = 1, 


and 

C 6 (t) = [f 7 (t) - P 6 (t>] • T ( t ) = 1 U(t) • u(t) = 1 


(38) 


which satisfies uniformity for 5 £ r £ 6 and yields a Polar system from 
the circular arc 


?(5,t) = 1 [^ 5 (t) + ? g (t)] = (j+ ft) u(t) , (39) 

up to the circular arc 

?(6,t) = fyt) = (f + /2) u(t). (40) 

The entire collection of bounding and intermediate surfaces are depicted 
in Fig. 8. 


96 








in the table. For a given mesh point, the interval k £ r < k + 1 
containing it determines the index k for G^, and G^-j respectively 
in Eqs. 22-24. Due to the uniform selection of partition points r^, a 
repetitive pattern in the evaluations can be observed and is 

indicative of translated versions of the same function. When 9 uniformly 
spaced mesh points are chosen for 0 £ t £ 1, and when the multi -surface 
transformation of Eq. 9 is evaluated for the 21 * 9 mesh, we obtain the 
coordinate mesh which is displayed in Fig. 9. From uniformity and 
Table 1, the first 8 and the last 5 mesh points in r are seen to be 
uniformly distributed, and the mesh is respectively purely Cartesian and 
purely Polar for those points. To illustrate the computational aspect, 
we shall explicitly evaluate the transformation at the point with curvi- 
linear variables r = 4.5, t = 0. At t = 0, we have 


3 ( 0 ) = ( 2 ( 0 ) - 1 , 0 ) = (- 1 , 0 ), 

and (41) 

u(0) = (-cos j, sin j) = (- -1, -~). 

q 4 /2 /2 

For r = 4.5, we are at the 15th mesh index in Table 1 where we read 
across to note that we are in the 4th interval (4 £ 4.5 £ 5) with 
G^(4.5) = .88 and G^(4.5) = .13. By substitution into the transformation 
(Eq. 25 for k = 4) we obtain 


S 
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^(4.5,0) = ? 4 (0) + G 4 (4.5)[? 5 (0) - ? 4 (0)] + G 5 (4.5)[? 6 (0) - ^(0)] 

= 5(0) + G 4 (4.5)[v^ u(0) - 5(0)] + G 5 (4.5) u(0) 

= (-1,0) + .88[/2(- \ 4) - (-1.0)] + . 13(- 4» 4) 

Jl /2 /2 /2 

= (-1,0) + .88(0,1) + (-.09,09) (42) 

= (-1 + 0 -.09, 0 + .88 + .09) 

= (-1.09, .97). 

ffe In continuation with local methods, the case with nonuniform parti- 

tions for the piecewise linear functions is given in Eiseman [5]. In 
addition, local interpolants with a higher level of smoothness (derivative 
continuity) can be used and are developed in Eiseman [7]. With the local 
controls over the transverse coordinate curves which connect two bounding 
surfaces, lateral bounding surfaces can also be approximately fit. A 
precise fit of the lateral boundaries can be obtained with blending 
functions which were used by Gordon and Hall [8] to create a global method. 
Further applications of blending functions will be presented at this work- 
shop by Ericksson [9], by Forsey, Edwards, and Carr [10] and by Anderson and 
Spradley [11]. 
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Algebraic Mesh Generation - Three Dimensions 


An algebraic approach to mesh generation in three dimensions results 
in algebraic functions that relate a computational domain to a physical 
domain. If the computational domain is defined by the three variables r, 
£> and z, on the unit cube 

0 £ r < 1 , 

0 < e < 1 , (43) 

0 < c < 1, 

then the physical domain in Cartesian (x,y,z) coordinates is given by the 
transformation ^(r,£,c) = (x,y,z) where 

x = x(r,E,c) , 


y = y(r,£»C) , (44) 

z = z(r,C»c). 

When Eq. 44 is nonsingular it has an inverse transformation denoted by 


r = r(x,y,z) , 


C = C(x,y,z) , ( 45 ) 

4 = c(x,y ,z) . 


A uniform mesh is defined on the computational domain by constants Ar, 
AC, A£ (Fig. 10). This mesh maps using Eq. 44 to a corresponding mesh 



Fig. 10 - Computational domain. 


in the physical domain which is not necessarily uniform. A simple example 
for Eq. 44 is given by 



( 46 ) 


or 

€ = 4 £n{l + (e kl - 1) f-}, 

K 1 \ 

r = in{l + (e 2 - 1)(^ - tan 6 y ) 1-}, 

? - j- tnfl + (e 3 - l)(f - tan e ) 4—} , 
k 3 x z l L 

where k-j , 1^, k^, 0^, e^, X^, Y^, and are constants. For 

Cq<_C£ 1» 0 £ r £ 1, 0<_^<1 and = Z^ the uniform computational 

domain maps into a frustrum of a paramid (Fig. 11) and the mesh is 
concentrated in the physical domain according to the magnitudes and 
signs of k-j , k£, and k^. 

Equation 44 must satisfy the constraints outlined in the 
introduction and which vary from problem to problem. For many mesh 
generation problems, the constraints reduce to having the boundaries 
in the computational domain map to boundaries in the physical 




Fig. 11 - Physical domain for Eqs. 46-47 < 


domain and concentrating the mesh in specified regions of the physical 
domain. The polynomial N-surface transformations (Eq. 6-18) are global 
algebraic mesh generation techniques which satisfy the basic boundary 
constraints and result in polynomial functions of degree (N - 1) 
with respect to one of the independent variables. For a small N the 
polynomials are particularly simple. If the surface coordinates are 
t = the transformation t^(r,t) = (x(r,t), y(r,t), z(r,t)) 



Fig. 12 - Boundary mapping 


s defined such that 


X D = x(0,F,»t) = X,(C.C)» 
B 1 

Y d = y(0,^jC) = Yj(£»£)» 
B 1 

Zn = Z(0,C»^) = Z-j(CjC) ) 
B 1 


Xn = X(lj5>C) 


x 2 (£ »£) » 


Yn = yO»£»C) = 
b 2 

Z B = z(l,5,e) = Z 2 (^,C), 



where = (X R , Y R , Z K ) is one boundary and "F* w (5,0 = 


B-j 5 'B 


'1 ’ B 1 

(Xg^, Yg , Zg ) is the other boundary in the physical domain (Fig. 12). 

The polynomial 4-surface transformation (Eq. 12 and 17) allows a 

constraint to be placed on the mesh in addition to that of fitting the 

boundaries. This constraint occurs when the physical mesh is required 

to be orthogonal at the boundaries. Since the derivatives ( 0 , c * c ) , 

~ (l,5,rj» etc. can be computed from the cross product of the tangential 
° r ’ dX-, dX, dY, dY, 

derivatives (£ ,c) , (£,c) , •^- L (C,c) , etc. , we have 


N' 


§£ U-l,f,c) t + fjr (fc-l,e,c) J + || U-1.C.0 $ = 


3r 


-* 

1 


J 


It 


5 




3Y £ 8Z £ 
?r^) 3 t ( c *° 


9X 0 3Y. 3Z. 

ar(^) ar(«»«) ir^) 


£ = 1,2 


(49) 


where T, (j, and it are unit vectors and K is the magnitude of the 
normal vector, the choice of which can be used to apply controls developed 
in Eiseman [1] for the shape of coordinate curves in r and for the 
distribution of constant r-surfaces. Applying this procedure will force 
the mesh to be orthogonal at the boundaries but not necessarily anywhere else 
A globally uniform computational mesh (for linear Sp in Fig. 6) can 
be mapped onto a physical mesh with the polynomial N-surface transformations 
given in Eqs. 6, 9, and 12. Concentration of mesh points in the r 
direction is accomplished by choosing a function r(r) such that 



Weigel [3] used the function 


r 



e 


k 


-1 


0 < r < 1 , 


(50) 


to contract the physical grid toward one boundary or the other. The number k 
is a free parameter whose magnitude dictates the degree of contraction. When r 
is replaced by r(r) in Eq. 6, the contractive function becomes embedded 
in the linear polynomial part of Eq. 6, which results in 
x = X 2 U,i;)r + X^S.eHl “ r). 
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y = Y 2 (£,<;)r + Y 1 (C,c)(l - r) , 
z = Z 2 (C»0' r ' ^ Z-|(^,^)(l " r) s 


(51) 


where a uniform partition of 0 £ r <_ 1 yields a desired nonuniform partition 
of 0 <_ r £ 1 that, in turn, proportionately partitions the linear segments of 
the transformation. 

The example previously presented can be derived with this approach 

where 


x(£) = X^C.C) = X 2 (5 ,c) = ?X L , 

= a L tan 0 y 

Y 2 U.5) = Y 2 (C) = a L (tan 9 y + YJ 


z(£,0 = Z-|(5»^) = Z 2 (£,<;) = £X L (tan 6 Z + 1^)1 + £X L tan 0 2 (1 - £) 

C = — c . (52) 


k 3^ 

e J -1 
k 3 

e J -1 


'1 


and 


x = £X l , 

y = Y 2 (Or + Y^OO - r), 

z = z(£,0> 
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where 



A fundamental constraint of this approach is the representation of 
the boundaries. The boundaries can be represented as analytical functions 
or approximate functions based on discrete data from the boundaries. 

In either case the representation must be in a form where parametric 
variables which can be normalized to the unit interval are the independent 
variables. If the parametric independent variables are chosen to be 
s and t, then for the two boundaries 

x^c.e) -► x-j ( s,t) , 

Y^.C) - Y-j (s ,t) , 

Z^c.c) - Z-j (s ,t) , 


X^^jC) "*■ X 2 (s ,t) , 
Y 2 (C,C) -*■ Y 2 (s,t), 
Z 2 (?»?) Z 2 (s ,t) , 


® ^ ^ ’ s min — s — s max’ 


oicii, t min 1 1 < t max . 


( 53 ) 
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The choice of parametric variables can vary from problem to problem. 


A relationship between (£>£) and (s,t) is 


s = C(s 
t = c(t 


max 


max 


s min^ + s min’ 

t . ) + t • . 
min' min 


(54) 


This is a linear relation which maps the unit interval onto the parametric 
variables. Contraction of the physical grid at the boundaries is 
accomplished in the same manner as for the internal grid distribution. 


I = C(0, > 0 , 

l = ^ > o, 

vJL, 

0 £ \ <_ 1 , 0 < £ < 1 , 
0 < \ £ 1 , 0 < £ £ 1 . 


(55) 


A pproximate Boundary-Fitted Coordinate Systems Using Tension Sp l ine 
Functions 

It is often the case that boundaries in a physical domain are 
described by discrete sets of points. The boundaries may be open or 
closed. An approximate boundary- fi tted coordinate system can be obtained 
using the technique described and a tension spline function interpolation 
to the discrete data defining the boundaries. Tension splines are chosen 
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because standard cubic splines and other higher order approximation 
techniques often result in wiggles in the approximation. Wiggles on a 
boundary using the technique propagate into the interior grid. The tension 
parameter embedded in the tension spline approximation to the curve 
allows control of the "curvedness" of the approximation. A very large 
magnitude of the tension parameter corresponds to a linear approximation, 
whereas a very small value corresponds to cubic splines. Tension splines 
can be applied in two and three dimensions. An example is presented 
here that is applicable to a two-dimensional mesh. 

Using the tension spline technique, a point set on boundary one is 
defined by {x^,y^}^~^ and on boundary two by {x^. ,yj } Approximate 
arc length is used as a parametric independent variable. The approximate 
arc length is: 


■ r ( \ + i 

0 

o 1/2 

- x i> + (j'i+i 

- y^ri 

■ r ( Vi 


O 1/2 

-v + <* J+ 1 

- V J 


i = 1 . . .n 


j = 1 . . .m 


S 1 


= 0 


(56) 


0 < s • < s 

— i - n 

0 < S- < S . 

— j — m 


no 


mil 


From the computational coordinate system the unit interval (0 < £, <_ 1) 
must be mapped onto each boundary; that is: 


s = s(0» 

0 <. £ <. 1 » 

0 <_ s <_ s n , 

0 < s < s m . 
— — m 


This is accomplished by letting 


( 57 ) 


s = £s n on boundary one and 
s = £s m on boundary two. 


(58) 


The tension spline functions are piecewise continuous hyperbolic 
functions on each boundary such that 


x = g"(s £ ) 


sinh[o(s £+1 - s)] 
o sinh[o(s £+1 - s £ )] 


+ 


+ 




sinh[o(s - s )] 


a sinh[a(s £+1 - s £ )] 


x * - 


g"(s t ) 


' S &+1 ~ s 
S £+l ‘ S £ , 



s - h 

Vl ' 2 

0 

\ S £+ 1 " S £ j 


(59a) 


m 



sinh[o{s , - s)] 
» - h "< s t > — 


cr sinh[a(s £+1 - s £ )] 


+ h"(s oxl ) 


sinh[a(s - s^)] 


£+ 1 ' 2 . 


a sinh[o(s £+1 - s £ )] 


y £ ' 




h “ (s £ ) 


S £+l - S 


S £+l S t 


h "<W 


S - 5, 


1S £+1 S £ 


V* 


(59b) 


£ = i on boundary one, 

£ = j on boundary two, 
a = tension parameter. 

The unknowns in these equations are g“ (s £ ) and h" ( s ^ ) which are second 

c — [ 

derivatives at the data points { 1 and are stained through enforce- 
ment of the continuity of the first derivatives at the data points, and 
the specification of two end conditions. A tri diagonal system of linear 
equations results for each set of unknowns. The solutions of the tridiagonal 
systems yield 9"(s„) and h"(s 0 ). 

e kr - 1 

A cubic polynomial and the contracting function r = e — r- - — provide 

e K - 1 

the relationship between the computational domain and physical domain. 


The deri vati ves 


dX, 

dn 


and 


dY < 

dn 


are: 


dn * ds ’ 

* dx * 

dn ds ’ 


(60) 
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' By defining a grid with constants A£ and Ar in the computational 

domain a corresponding grid is explicitly defined in the physical domain. 
§. An example of an airfoil grid is presented (Fig. 14). Data points on 

- each boundary, magnitude of the normal derivative, and contract parameter 
values define the grids. Boundary data for the airfoil is shown in 
Fig. 15. Also, for closed boundaries such as the airfoil, periodic 
conditions are applied in the spline functions. 



Fig. 14 - Mesh for Karman Trefftz Fig. 15 - Data definition for 

airfoil using splines under boundaries of Karman Trefftz 

tension. airfoil domain. 


Complex Three-Dimensional Mesh Generation 


The development of three-dimensional meshes where mesh lines are free 


to move in all three coordinate directions is extremely difficult. The 
reader is directed to reference 9 for examples of such unconstrained three- 
dimensional meshes. An expedient approach for complex three-dimensional 
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geometries is to attempt to simplify the problem by rendering the three- 
dimensional geometry quasi two-dimensional. This is the essence of the 
frustrum pyramid mesh previously presented. Also when there is axis- 
symmetry in a problem, two-dimensional rendering of the geometry is 
possible. This is demonstrated with the spike-nosed configuration shown 
in Fig. 16. Figure 17 shows an algebraic generated mesh for one plane 



Fig. 16 - Surface for a spike-nosed body. 


of the geometry. The mesh is made three dimensional by rotating the mesh 
in movements about the axis of symmetry. 

For aircraft surfaces the problem is more difficult. There are, 
however, several good software packages ([12] is a good example) for 
surface definition of aircraft aerospacecraft geometries. In [12] 




Fig. 17 - Mesh for spike-nosed body. 


the Harris geometry is used to establish coordinates for Coon's surface 
definition of a configuration. Spline functions are used to compute 
the derivatives for the Coon's surface definition. Figure 18 shows the 
input description for a wing-fuselage configuration where the wing and 
fuselage are defined separately. Figure 19 shows an enriched definition 
of the configuration using the Coon's surfaces. A part of the available 
software described in [12] is the ability to compute the intersection 
of an arbitrarily defined plane and the surface definition. For the 
configuration shown in Fig. 19 planes perpendicular to the x axis and 
at a constant increment in the x direction are shown in Fig. 20. 

If an outside boundary is defined the corresponding intersection with the 
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Fig. 20 - Planar intersections with a wing-fuselage configuration. 


planes is computable. With an inside and outside boundary the techniques 
previously described can be employed in two dimensions. Complexities 
of multi connected regions is introduced but this can be attacked with 
branch cuts (Fig. 21). 

ORIGINAL PAGE IS 
O v pO'”’ 



Fig. 21 - Branch cuts for multi -connected sections. 


Another attack on complex three dimensional problems is to define 
several computational domains (Fig. 22) with mutual boundaries. The mapping 



Fig. 22 - Multiple computational domains with mutual boundaries. 


function transforms the computational domains to parts of the physical 
domain with mutual intersections. 

Conclusions 

Algebraic methods provide precise controls for mesh generation. 
Methodologies for mesh construction can be based on a parameterized 
description of surfaces which consist of bounding surfaces and intermediate 
control surfaces. The surface locations determined by the respective 
surface parameterizations determine the nature of the transverse coordinate 
curves which connect the bounding surfaces. Relative to uniform conditions, 
precise control over the mesh placement in the physical domain can be 


accomplished by embedding distribution control functions in the surface 
parameterizations or in the transverse direction. Complex bounding 
topologies, especially in three dimensions, cause mesh construction 
difficulties. It is proposed that whenever feasible, the complex topol- 
ogy be simplified such as rendering the geometry quasi two-dimensional. As 
a final remark, precise controls are one of the major advantages of algebra 
methods: they give the capability to prescribe specific desirable and 

helpful mesh formations. 
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